program myprobit, sortpreserve
	version 8.1
	if replay() {
		if (`"`e(cmd)'"' != "myprobit") error 301
		Replay `0'
	}
	else	Estimate `0'
end

program Estimate, eclass
	syntax varlist [if] [in]		///
		[fweight pweight] [,		///
		noLOg noCONStant		/// -ml model- options
		Robust CLuster(varname)		///
		noLRTEST svy			///
		OFFset(varname numeric)		///
		EXPosure(varname numeric)	///
		Level(integer `c(level)')	/// -Replay- option
		*				/// -mlopts/svyopts- options
	]

	// mark the estimation sample
	marksample touse

	// check syntax
	if "`svy'" != "" {
		svymarkout `touse'
		local wvar `r(weight)'		// remember weight var
		svyopts modopts diopts options , `options'
		if "`s(subpop)'" != "" {
			local subpop `s(subpop)'
			local subopt subpop(`subpop')
		}
		local meff `s(meff)'
	}
	local diopts `diopts' level(`level')
	mlopts mlopts, `options'
	local cns `s(constraints)'
	gettoken lhs rhs : varlist
	if "`cluster'" != "" {
		local clopt cluster(`cluster')
	}
	if "`weight'" != "" {
		tempvar wvar
		quietly gen double `wvar' `exp'
		local wgt "[`weight'=`wvar']"
	}
	if "`wvar'" != "" {
		// for initial value calculations
		local awgt "[aw=`wvar']"
	}
	if `"`robust'`cluster'`cns'"' != "" | "`weight'" == "pweight" {
		local lrtest nolrtest
	}
	if "`log'" != "" {
		local qui quietly
	}
	if "`offset'" != "" {
		local offopt "offset(`offset')"
	}
	if "`exposure'" != "" {
		local expopt "exposure(`exposure')"
	}

	// markout missing values from the estimation sample
	markout `touse' `wvar' `offset' `exposure' `subpop'
	markout `touse' `cluster', strok
	if "`subpop'" != "" {
		// restrict initial value calculations to
		// observations within the subpopulation
		local subtouse `touse' & `subpop' != 0
	}
	else	local subtouse `touse'

	if "`meff'" != "" {
		if "`constant'" == "" {
			// initial value for the misspecified model
			sum `lhs' if `subtouse', mean
			local b0 = invnorm(r(mean))
			local initopt init(_cons=`b0') search(quietly)
		}

		`qui' di as txt _n "Fitting misspecified model:"
		// fit the misspecified model
		ml model d2 myprobit_d2				///
			(xb: `lhs' = `rhs',			///
				`constant' `offopt' `expopt'	///
			)					///
			if `subtouse',				///
			`log' `mlopts' `initopt'		///
			missing maximize
		// covariance matrix from the misspecified model fit
		tempname Vmsp
		matrix `Vmsp' = e(V)
	}

	if "`constant'" == "" {
		// initial value
		sum `lhs' `awgt' if `subtouse', mean
		local b0 = invnorm(r(mean))
		local n = r(N)
		local initopt init(_cons=`b0') search(quietly)

		if "`lrtest'" == "" {
			quietly count if `subtouse' & `lhs' != 0
			local n1 = r(N)
			local lf0 = `n1'*ln(norm(`b0'))		///
				+(`n'-`n1')*ln(norm(-`b0'))
			local initopt `initopt' lf0(1 `lf0')
		}
	}

	// fit the full model
	`qui' di as txt _n "Fitting full model:"
	ml model d2 myprobit_d2				///
		(xb: `lhs' = `rhs',			///
			`constant' `offopt' `expopt'	///
		)					///
		`wgt' if `touse',			///
		`log' `mlopts' `robust' `clopt'		///
		`svy' `modopts' `initopt'		///
		missing maximize

	if "`Vmsp'" != "" {
		// compute misspecification effects and save them to e()
		_svy_mkmeff `Vmsp'
	}
	// save a title for -Replay- and the name of this command
	ereturn local title "My probit estimates"
	ereturn local cmd myprobit

	Replay , `diopts'
end

program Replay
	syntax [, Level(integer `c(level)') * ]
	// use svyopts to get display options
	svyopts noopts diopts, `options'
	// anything in `noopts' will cause an error
	mlopts noopts, `noopts'
	ml display , level(`level') `diopts'
end
